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Recent concrete proposals suggest it is possible to engineer a two-dimensional bulk phase support¬ 
ing non-Abelian Fibonacci anyons out of Abelian fractional quantum Hall systems. The low-energy 
degrees of freedom of such setups can be modeled as Z 3 parafermions “hopping” on a two-dimensional 
lattice. We use the density matrix renormalization group to study a model of this type interpolating 
between the decoupled-chain, triangular-lattice, and square-lattice limits. The results show clear 
evidence of the Fibonacci phase over a wide region of the phase diagram, most notably including 
the isotropic triangular lattice point. We also study the broader phase diagram of this model and 
show that elsewhere it supports an Abelian state with semionic excitations. 


I. INTRODUCTION 


The experimental realization of non-Abelian anyons— 
emergent particles with highly exotic exchange 
statistics—is anticipated to have widespread impli¬ 
cations. Apart from demonstrating a profound new 
facet of quantum mechanics, many applications await 
the completion of this ongoing quest. Examples include 
tests of Bell’s inequalities,® robust quantum memory, 
novel low-temperature circuit elements,® and most 
importantly intrinsically fault-tolerant ‘topological’ 
quantum computing.®® The last of these relies on the 
remarkable fact that adiabatic exchange (braiding) of 
non-Abelian anyons enacts a unitary rotation within 
the space of locally indistinguishable ground states 
generated by the anyons. Storage and manipulation of 
qubits encoded in these ground states thus takes place 
non-locally, so that the quantum information is securely 
‘hidden’ from local environmental perturbations. 

While most of these applications can be carried out 
with any non-Abelian anyon type, topological quan¬ 
tum computation carries more stringent demands. Con¬ 
sider Ising anyons—or defects that bind Majorana zero 
modes—which likely constitute the most experimen¬ 
tally accessible non-Abelian anyon. Numerous realis¬ 
tic platforms now exist for realizing I sing anyons, most 
prominently in quantum Hall systems^ and topologi¬ 
cal superconductorsand indeed tantalizing experi¬ 
mental evidence of these particles has accumulated in 
both settings.ll^lf^ Braiding Ising anyons, however, only 
amounts to 90° qubit rotations on the Bloch sphere. Per¬ 
forming the arbitrary qubit rotations necessary for uni¬ 
versal computing with Ising anyons requires additional 
operations that are not topologically protected. 

This shortcoming strongly motivates the pursuit 
of other types of non-Abelian anyons with ‘denser’ 
braid statistics. One appealing alternative class are 
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FIG. 1. Model of Z 3 parafermions with interactions form¬ 
ing an anisotropic triangular lattice with intrachain couplings 
ts (horizontal) and interchain couplings ti (vertical) and t 2 
(slanted). 


defects binding parafermionic zero modes^ which 
comprise natural Majorana generalizations. Al¬ 
though such defects require a strongly interacting 
host system (unlike Ising anyons), many plausi¬ 
ble realizations have been suggested such as lat¬ 
tice quantum anomalous Hall states,!^ Abelian quan¬ 
tum Hall/superconductor het erostr ucturesmul¬ 
tilayer quantum Hall systemsand coupled-wire 
arrays.l^^^®^ Parafermionic zero modes produce a larger 
ground-state degeneracy than Majorana modes and thus 
enable a denser set of qubit rotations through braid¬ 
ing. While providing some advantages for quantum 
computation,!^ their braid statistics nevertheless remains 
non- universal.!^ 

Fortunately one can leverage setups supporting 
parafermionic zero modes to generate non-Abelian 
anyons allowing bona fide universal topological quantum 
computation. To illustrate how this is possible, imagine 
nucleating a two-dimensional (2D) array of parafermionic 
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zero modes in, say, a quantum Hall/superconductor hy¬ 
brid structure. The collection of zero modes encodes 
a macroscopic ground-state degeneracy, similar to a 
partially filled Landau level. Hybridizing parafermions 
among the sites can lift this degeneracy and produce new, 
possibly very exotic 2D phases. Resolving the precise 
state selected by the coupled parafermion modes, how¬ 
ever, poses a highly nontrivial task since the system can 
not be described by free particles (much like the Landau 
level problem with Coulomb interactions). 

Inspired by the important work of Teo and Kane,S^ 
Ref. 03] nevertheless identified an analytically tractable 
limit where strongly anisotropic, weakly coupled 
parafermion chains could be shown to enter a topolog¬ 
ically ordered ‘Fibonacci phase’. This phase represents 
a cousin of the so-called Z 3 Read-Rezayi non-Abelian 
quantum Hall state,^^ yet is built from well-understood 
Abelian states of matter. (The analogy closely parallels 
the relation between a p + ip superconductor and the 
Moore-Read state.^^^^ Like its quantum Hall cousin, the 
Fibonacci phase supports Fibonacci anyons, which are 
the holy grail for topological quantum computation since 
they allow one to approximate arbitrary unit ary gat es to 
any desired accuracy solely through br aiding .E®46| 

While this progress is encouraging, the stability of the 
Fibonacci phase away from the solvable quasi-ID limit— 
as well as the broader phase diagram for coupled 2D 
parafermion arrays—remains largely unknown. Obtain¬ 
ing a more thorough and quantitative understanding of 
these physical systems poses a pressing issue given the 
importance of finding experimentally accessible realiza¬ 
tions of Fibonacci anyons for quantum computing. In 
this paper we take a major step in this direction by per¬ 
forming extensive density matrix renormalization group 
(DMRG)I^^®^ simulations of lattice parafermions with 
couplings sketched in Fig. A virtue of the model we 
study is that it interpolates between various physically 
interesting regimes, including decoupled chains and the 
isotropic triangular- and square-lattice limits as special 
cases. 

DMRG naturally complements the weakly-coupled 
chain approach to realizing the Fibonacci phase. A 
key observation from analytics is that the correla¬ 
tion length for the Fibonacci phase becomes arbitrarily 
small with weakly coupled chains, albeit at the expense 
of extended ‘x’ correlations. This tradeoff is highly fa¬ 
vorable for numerics. The short y correlation length al¬ 
lows us to well-approximate the 2D limit of interest us¬ 
ing systems composed of relatively few chains. Moreover, 
with sufficient effort DMRG can handle very long chains 
in the x direction even when each is close to criticality, 
facilitating direct comparisons to analytical predictions. 
Near the decoupled-chain limit we can unambiguously 
pinpoint the Fibonacci phase in our simulations, then 
systematically track how it evolves as we gradually en¬ 
hance the interchain coupling—all the while keeping the 
y correlation length manageable. The latter feature is 
important since in practice accurate DMRG results for 


typical 2D models on cylinders can only be obtained for 
circumferences on the order of ten sites along y. Finally, 
we note that moving away from the decoupled-chain limit 
lessens the overall entanglement, mitigating the cost of 
the additional sites needed along the y direction and al¬ 
lowing DMRG to even address the behavior of isotropic 
2D systems. 

In this way we numerically show that the Fibonacci 
phase extends across a very wide swath of parameter 
space in our model, persisting from weakly coupled chains 
all the way to the isotropic triangular lattice point and 
beyond. Evidently the quasi-lD limit pursued earlier is 
by no means necessary, but merely provides a convenient 
entry way into the relevant physics. The broad stability 
window for Fibonacci anyons that we identify is one of the 
main punchlines of this paper. We also present evidence 
for a second topologically ordered state that is Abelian 
and supports a semion as its only nontrivial quasiparti¬ 
cle. This phase naturally emerges from weakly coupled 
chains upon swapping the sign for the intrachain cou¬ 
pling constant, and in our simulations also appears quite 
robust even away from this limit. In fact the Fibonacci 
and semion phases comprise the only two states that ap¬ 
peared throughout the broad (though not completely ex¬ 
haustive) parameter space that we explored numerically. 
For a summary of our results see the octahedral phase 
diagram presented as a cutout template in Fig. 

We organize the remainder of the paper as follows. 
Section [H] motivates the lattice model we study from 
the viewpoint of physical quantum Hall-based architec¬ 
tures. Section [Hi] analyzes the model, both analytically 
and numerically, in the two-chain limit—which we argue 
already contains precursors of Fibonacci physics. The 
multi-chain setup of greatest interest is tackled in Sec. HYi 
We present a wealth of numerical evidence indicating the 
robustness of the Fibonacci phase; analytically capture 
the semion state; and explore the broader phase diagram 
of the model. Section m discusses future directions and 
highlights the connection between our study and recent 
related works. Four appendices contain additional details 
of our model and simulation methods. 


II. Z3 PARAFERMION LATTICE MODEL 

To motivate the lattice model we will study, let us re¬ 
view how a non-Abelian phase supporting bulk Fibonacci 
anyons can arise from coupling parafermionic zero modes 
nucleated in a fractional quantum Hall fluid. As a 
primer consider some arbitrary Abelian quantum Hall 
state sliced into two adjacent halves as shown in Fig.j^a). 
The cut produces a new set of gapless counterpropagat- 
ing edge states opposite the trench. One can always fully 
gap these modes—effectively resewing the two halves—in 
more than one physically distinct way. The most natu¬ 
ral mechanism involves backscattering electrons across 
the trench to simply recover the original uninterrupted 
quantum Hall state as Fig. |^b) illustrates. Filling the 
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FIG. 2. (a) A trench dividing two regions of an Abelian quan¬ 
tum Hall phase supports gapless modes, (b) Gapping these 
modes via electron backscattering recovers the original bulk 
phase, (c) Introducing either Gooper pairing, or ‘crossed’ 
tunneling in a bilayer system, provides an inequivalent way of 
gapping the trench, (d) Domain walls between incompatibly 
gapped regions support protected zero modes. 


trench with a superconductor [blue region in Fig. ic)] 
provides a second, intuitively quite different method: the 
edge modes can then gap out by Cooper pai ring elec¬ 
trons from opposite sides of the trench.^^^J^ Some se¬ 
tups can support alternative charge-conserving gapp ing 
mechanisms (aside from trivial backscattering) Por 
instance, in quantum Hall bilayers formed out of Laugh- 
lin states at filling v = 1 /m, one can gap the trench by 
‘crossed’ tunneling whereby electrons hop between the 
top layer on one side and the bottom layer on the other, 
sewing the halves in yet a different manner. Domain walls 
[Fig. [^d)] separating regions gapped in such incompati¬ 
ble ways bind protected zero-energy modes that form the 
basic building blocks in our lattice model. 

We are specifically interested in setups that support 
Z 3 parafermionic zero modes, which provide the min¬ 
imal extension of the more familiar Major ana opera¬ 
tors. Concrete realizations include (i) domain walls be¬ 
tween pairing- and tunneling-gapped trenches in spin- 
unpolarized z/ = 2/3 quantum Hall fluidand (ii) do¬ 
main walls between trivial and ‘crossed’ electron tunnel¬ 
ing in = 1/3 bilayers.l^ Consider for the moment a one¬ 
dimensional array of well-separated Z 3 zero modes as in 
the geometry displayed in Fig. These modes encode a 
ground-state degeneracy that can be understood as fol¬ 
lows. In the z^ = 2/3 realization each superconducting- 
gapped region can accommodate charge 2e/3 without en¬ 
ergy penalty, while in the bilayer example each ‘crossed’ 
region can similarly acquire an e/3 dipole for free. 

Parafermionic zero-mode operators cycle the system 
through the corresponding degenerate ground-state man¬ 
ifold. Importantly, charges can be added to the trench 
in two physically distinct ways—and hence there exists 
two inequivalent representations of the parafermion oper¬ 
ators. We denote these two parafermion representations 



FIG. 3. One-dimensional chain of parafermionic zero modes 
arising from inequivalently gapped regions of a trench. Hy¬ 
bridization among nearest-neighbors in the chain is described 
by the Hamiltonian in Eq. 


by OiR ^2 oll ,3 with j the domain-wall site index. As a 
specific example, in the z^ = 2/3 setup and al¬ 
ter the adjacent superconducting region by adding charge 
2e/3 to the upper and lower trench edges, respectively. 
This distinction is meaningful since fractional charge can 
not migrate across the trench. For further discussion on 
this important point see Ref. [43] as well as Appendix 
Both representations fulfill the conditions 

= 1; A,j = (1) 

for A = R or L; moreover, they exhibit ‘anyonic’ com¬ 
mutation relations 


aRj^Rj’ = e 
OlL,jOlL,j' = e 




( 2 ) 


inherited from the quantum Hall edge fields from which 
they derive. 

Suppose now that the domains in Fig.j^are sufficiently 
narrow that nearest-neighbor parafermionic modes hy¬ 
bridize appreciably. One can describe the hybridization 
with a Hamiltonian 


where uj = n sums over unit cells, and the J and 

h terms couple alternating pairs; see Fig. Physically, 
these terms reflect tunneling of fractional charges from 
one domain wall to the next via a path above the trench. 
(Note that an equivalent form in terms of aR^j is also 
possible; this form would correspond to tunneling paths 
below the trench.) 

As reviewed in Appendix the single-chain Hamilto¬ 
nian in Eq. (§ maps to the three-state quantum Potts 
model under a non-local ‘Fradkin-Kadanoff’ transforma¬ 
tion a kin t o the Jordan-Wigner mapping in the Ising 
mo del.Many properties immediately follow from this 
identification. Like its Potts analogue, with J, /z > 0 
the Hamiltonian admits two distinct phases. For h > J 
the trench enters a trivial state where all parafermion 
modes dimerize and gap out pairwise with their neigh¬ 
bors. On the other hand, for J > h the parafermion 
operators dimerize in a shifted pattern leaving ‘unpaired’ 
modes at the ends of the trench (assuming open boundary 
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FIG. 4. A series of one-dimensional parafermion chains inter¬ 
act by tunneling fractional charge through the quantum Hall 
fluid [Eq. @i- For the present work, we consider anisotropic 
next-neighbor tunneling processes for parafermions arranged 
in a triangular lattice, allowing us to study (among other lim¬ 
its) the isotropic triangular lattice {ti = t 2 = ^ 3 ), the isotropic 
square lattice (ti = ts, t 2 = 0), and weakly coupled chains 
{ts > ti, ^ 2 ). 


conditions).!^ Here the system exhibits a protected three¬ 
fold ground-state degeneracy similar to the topological 
phase of a Kitaev chain.!^ These two phases are separated 
by a self-dual critical point at = J described by a non- 
chiral Zs-parafermion c onfor mal field theory (CFT) with 
central charge c = 4/51111111 As our nomenclature hints 
the right- and left-moving parafermion fields at this crit¬ 
ical point are closely related to the lattice parafermion 
operators <yR/Lj^ 

Now consider a series of N such trenches at vertical 
positions y, arranged such that parafermion modes re¬ 
side at the sites of a 2D lattice (Fig.[^. If trenches y and 
^ +1 are separated by a finite distance, there will also be 
fractional-charge tunneling processes through the inter¬ 
vening quantum Hall fluid. Importantly, these processes 
can only couple the deft’ parafermion representation in 
trench y to the ‘right’ representation in ^ + 1. (The corre¬ 
sponding process with right and left interchanged would 
pass fractional charge through regions where only whole 
electrons can pass.) A general interchain Hamiltonian 
thus reads 

N-l 

^- = -EE[ + 1) + H.c.j . 

?/=i 

(4) 

Because the bulk of each quantum Hall fluid is gapped, 
the interchain couplings tj-j> should decay rapidly as a 
function of separation between sites. Thus taking tj-j/ 
non-zero for only the first- or possibly second-nearest- 
neighbors in adjacent trenches constitutes a reasonable 
minimal model. 

Figure depicts the specific pattern of interactions 
that we consider. Parafermion modes bound to domain 
walls reside at the sites of a triangular lattice and couple 
via anisotropic next-neighbor tunneling processes. The 


Hamiltonian of this model reads 

H = Ho 

^o = -EE ^3 

y 3 

^- = -EE[ tiua^^^^{y)aRj{y + 1) '' ^ 

y 3 

Our intrachain Hamiltonian Hq is identical to Eq. § 
but fixed to J = h = t^ such that each horizontal 
parafermion chain remains tuned to its critical point if 
decoupled from the other chains. The interchain bonds 
then form the triangular lattice pattern of interest, with 
vertical couplings ti and slanted couplings ^2 in the 
skewed layout of Fig. Throughout we will assume 
purely real Importantly, the Hamiltonian’s phase 

diagram is sensitive to both the magnitude and signs for 
for these couplings, since the only gauge transformations 
we can make correspond to shifting the parafermion oper¬ 
ators by factors of uj. Notice that the model interpolates 
between the square lattice (^2 = 0 , ti = ts), isotropic 
triangular lattice (ti = t 2 = ts), and decoupled chain 
(ti = ^2 = 0 ) limits—allowing us to explore several im¬ 
portant cases in one framework. 

Despite the fact that Eq. (§ describes an inherently 
strongly interacting 2D setup, one can make controlled 
analytic progress in the weakly-coupled-chain limit with 
0 < ti = ^2 ^ ts. In this regime the problem reduces to 
parafermion CFTs in each chain that hybridize to yield 
a 2 D topologically ordered ‘Fibonacci phase’.(See Ap¬ 
pendix 1^ for details.) Quite generally, the Fibonacci 
phase exhibits the following key properties: 

{i) Within the parent quantum Hall fluid, spatial 
boundaries of the Fibonacci phase support ehiral Z 3 
parafermion modes with central charge c = 4/5. These 
residual gapless edge modes represent ‘half’ of the non- 
chiral CFT that occurs in each critical chain. 

(a) The bulk supports two quasiparticles types: trivial 
(fermionic/bosonic) excitations and Fibonacci anyons. 

(in) If the system is defined on a torus or an infinite 
cylinder, there are two degenerate ground states ll) and 
\e) —one associated with each quasiparticle type.l^ 

{iv) The entanglement entropy of ground state 
\n = l^s) on an infinite cylinder of large but finite 
circumference Ny scales as Sn = aNy — 7 ^. Here 
jn = \og{V/dn) represents the topological entanglement 
entropy expressed in terms of the anyon quantum dimen¬ 
sions = 1 , dg = (/9 and the total quantum dimension 
V = with ip = {1 Vb)/2 the golden ratio. 

(v) The entanglement spectrunP^l —which can be com¬ 
puted from the CFT!^^! —resembles the energy spectrum 
of the physical edge with chiral central charge 4/5. 
The entanglement spectrum of |1) consists of states de¬ 
scended from primaries (with scaling dimen¬ 

sions {0, |, |} respectively), while the spectrum of \s) 


«ki+i(y)«fl,i(y) + H-C. 
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consists of descendents of {s,a,a^} (with scaling dimen¬ 
sions {§5 respectively). The specific countin^^of 

states appears in Fig. [m 

Outside of the highly anisotropic, weakly-coupled- 
chain limit, the analytical methods used to establish the 
preceding results break down entirely. Addressing the 
broader phase diagram of the model—particularly the 
extent of the Fibonacci phase—is the central goal of this 
paper. For this we turn to density matrix renormaliza¬ 
tion group (DMRG) calculations of the ground states of 
Eq. (§, primarily on infinitely long cylinders using the 
infinite DMRG algorithm proposed in Ref. |58l To imple¬ 
ment the Hamiltonian in Eq. © we do not work directly 
with parafermion operators [Eqs. 0 ©]■ Instead, fol¬ 
lowing Appendix!^ we map the parafermion degrees of 
freedom to Z 3 clock variables, which provides a much 
more convenient (but formally equivalent) representation 
for numerics. 

It is worth noting some technical features of the model 
in Eq. 0 that facilitate our DMRG studies. Retain¬ 
ing all possible first- and second-neighbor interchain cou¬ 
plings in a square-lattice arrangement of the parafermion 
sites would lead to a one-dimensional Hamiltonian (as 
seen by the DMRG algorithm) with interactions up to 
a range dmax = where Ny is the circumference of 

the quasi-2D cylinder used for the simulation. In con¬ 
trast, the pattern of triangular lattice couplings we study 
here (Eig. gives a one-dimensional Hamiltonian with 
a maximum range scaling only as dmax = Ny. Addition¬ 
ally, as we discuss in the next section, the two-chain limit 
{Ny = 2 ) of the triangular model with open ^-boundary 
conditions is self-dual for a certain choice of parameters. 
Appealing to self-duality allows us to exactly determine 
the location of an important line of critical points that 
roughly represents a remnant of the Eibonacci phase com¬ 
pressed into a two-chain system. Einally, the properties 
of the bona fide 2 D Eibonacci phase enumerated above 
provide sharp numerical fingerprints that we can use to 
track the phase diagram in our multi-chain simulations 


{Ny > 2 ), which we describe in Sec. IV 


III. TWO-CHAIN LIMIT 


Before considering the model Eq. ^ on multi-leg cylin¬ 
ders, it is helpful to understand the limit of only two 
coupled chains with open boundary conditions in the y 



FIG. 5. Two-chain limit of the Hamiltonian in Eq. ([^, with 
open boundary conditions along the y direction. 


<•'111111 
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FIG. 6. Pairing of sites in the ground state for the two-chain 
ladder system in the limits of (a) dominant ti >0, which 
produces a trivial gapped phase, and (b) dominant t 2 > 0, 
which produces a topological gapped phase with a three-fold 
ground-state degeneracy. 


direction (see Fig. [^. As we will see this limit already 
contains precursors of the 2 D Fibonacci phase that we 
will uncover later on when studying wider geometries. 
Throughout this section we assume ts > 0 for simplicity 
but consider either sign of ti^ 2 - 

To understand the phase diagram of the two-chain sys¬ 
tem, first consider the limit where ti > 0 greatly exceeds 
both t 2 and ts. In the extreme case with ^2 = ts = 0, 
the ground state is found by pairing each parafermion 
site with the one directly above or below it, yielding the 
trivial product state illustrated in Fig.j^a). (Saying that 
two parafermions ai and aj ‘pair’ means that they form 
an eigenstate of ^(cj ato^j + H.c.) with maximal eigen¬ 
value 1.) Because this limit supports a unique ground 
state protected by a robust gap, this trivial phase per¬ 
sists upon restoring sufficiently small ^ 2,3 couplings. 

Gonsider next the limit with ^2 > 0 much larger than ti 
and ts. With ti = ^3 = 0 the ground state again arises by 
pairing sites, but now in the skewed pattern depicted in 
Fig.§b). Any finite ladder with open boundaries along 
the horizontal direction thus contains one unpaired site 
at each end—signifying a topologically nontrivial phase. 
These two decoupled end sites, taken together, form a 
degenerate three-level system. It follows that there are 
three degenerate ground states distinguished by their ‘tri- 
ality’, defined as 


Q = , (6) 

n 

which admits eigenvalues 1, cj, or Restoring small, 
finite ti and ts only splits the ground-state degeneracy by 
an exponentially small amount in the system size, as the 
processes mixing the end states require tunneling across 
a macroscopic number of sites. (Note that this argument 
applies only to the ground states. The excited states, 
strictly three-fold degenerate for ti = ts = 0 , generally 
split by an amount decaying only as a power law in the 
system size.l^ 

Intuitively, one expects a ID phase transition between 
the gapped states sketched in Fig. when ti and ^2 be¬ 
come comparable. To be more quantitative, we invoke a 
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FIG. 7. Bulk gap of the two-chain system (Fig. as a func¬ 
tion of 0 with ti = sin^, t 2 = cos^, and ta = 1- The vertical 
dashed line marks the ti = t 2 (= l/\/2) critical point separat¬ 
ing topological and trivial gapped phases. For the topological 
phase, the first excited state is in the quantum number sector 
with triality 1; for the trivial phase, the first excited states 
are two-fold degenerate and lie in the u or triality sectors. 


weakly-coupled-chain analysis and find an interesting sce¬ 
nario. In the strictly decoupled-chain limit (ti^2 = 0 ), the 
low-energy properties of each chain can be described by 
a Z3 parafermion CFT at c = 4/5 with a pair of counter- 
propagating chiral fields . 1 ^ Turning on weak interchain 
couplings ti^2 generically gaps all of these fields, but when 
one fine-tunes ti = ^2 > 0 the top chain’s left-mover and 
the bottom chain’s right-mover remain gapless.^^ Thus 
along the line ti = ^2 (at least for 0 < ti,t2 ts) the 
two-chain system should be critical and described by a 
single non-chiral Z3 parafermion CFT in which the right- 
and left-movers have, in a sense, spatially separated along 
the vertical direction. We can, however, actually reach 
a much stronger conclusion. Appendix 0 shows that for 
ti = t2 the two-chain Hamiltonian maps to itself un¬ 
der duality followed by a time-reversal transformation. 
Hence there should exist a continuous phase transition 
line precisely along ti = ^2 > 0 for any ts > 0. 

We confirm all of the above predictions—and ex¬ 
tend them beyond the analytically tractable regimes— 
using DMRG. Unless otherwise stated we use the infinite 
DMRG algorithnJ^ to reach the thermodynamic limit 
in the horizontal direction while taking open boundary 
conditions along the vertical direction. For both gapped 
phases we find essentially exact results (truncation er¬ 
rors below 10“^^) by keeping only a few hundred states 
in DMRG. First, we check for the presence of the criti¬ 
cal line at ti = ^2 by computing bulk gaps using infinite 
boundary conditions.!^ As shown in Fig. the gap in¬ 
deed closes precisely at ti = ^2- (In the horizontal axis 
0 parameterizes the couplings through ti = sin^ and 
t2 = cosO with ts = 1.) 

To additionally extract the central charge along this 


FIG. 8. Entanglement entropy Sl{x) for a length-x subregion 
of a periodic two-chain system with length L = 60. Data 
were obtained at the critical point using parameters ti/ts = 
^ 2/^3 = 1- By fitting to the GFT prediction we extract a 
central charge c = 4/5, consistent with analytical arguments. 


critical line, we use DMRG to study periodic systems of 
size L along the horizontal direction and compute the en¬ 
tanglement entropy Sl{x) of a subregion of size x. For a 
CFT with central charge c the entanglement is predicted 
to scale a^ 


Sl{x) = 3 In 



( 7 ) 


Figure displays our simulation results. For ti/t^ = 
^2/^3 = 1 we find c = 4/5 to very high accuracy, con¬ 
firming the prediction that, along the critical ti = t2 
line, the coupled chains host a single pair of gapless Z3 
parafermion CFT edge modes. 

Figure shows the full ^3 > 0 phase diagram of the 
two-chain system computed with DMRG. We numeri¬ 
cally identify the trivial gapped phase occurring in the 
large ti > 0 region by searching for the presence of a 
finite gap and a unique ground state. To identify the 
region lying in the topological phase stabilized for large 
t2 > 0, we primarily examine the entanglement spectrum 
within DMRG. For the entire topological phase each en¬ 
tanglement “energy” exhibits a robust three-fold degen¬ 
eracy. At various points within this phase we also checked 
that finite systems with open boundary conditions pos¬ 
sess three degenerate ground states. 

In addition to the above two gapped phases, we also 
find a large region for either ti or ^2 negative where the 
system is gapless. Numerically computing the central 
charge using the same approach as for Fig.[^ we find this 
entire gapless region has central charge c = 1 . To under¬ 
stand this phase, consider the limit ti < 0 and ^2,^3 ^ 0. 
In this limit the antiferromagnetic ti coupling favors two 
of the three states on each rung, effectively projecting 
out the third. This allows a natural mapping to a spin 
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FIG. 9. Full phase diagram of the two-chain system in units 
where ta = 1. The trivial phase (green circles) surrounds the 
> 0 axis while the topological phase (blue triangles) sur¬ 
rounds the t 2 > 0 axis. The remainder of the phase diagram 
shown is gapless with central charge c = 1, which we verified 
numerically for points indicated with red squares. Transition 
lines separating the gapless state from the two gapped phases 
are approximate. The critical point along the line ti = t 2 is 
exact, however, and represents a precursor to the Fibonacci 
phase appearing in the 2D limit. 


1/2 chain. Perturbatively reintroducing t2 and ts gives a 
Hamiltonian of XXZ type. Close to the negative ti axis, 
the effective Z coupling is parametrically suppressed rel¬ 
ative to the XY coupling, leading to a gapless XY phase 
that is known to have c = 1. 

For our purposes the most interesting feature of the 
two-chain phase diagram is the c = 4/5 transition line 
occurring at ti = t2 where, remarkably, the low-energy 
right- and left-movers residing at opposite sides of the 
ladder do not hybridize even for large ti^2- This criti¬ 
cal point reflects a maximally ‘squeezed’ cousin of the 
stable 2D Fibonacci phase; in the latter, the right- and 
left-movers are separated by macroscopic distances and 
therefore coexist harmoniously without any fine-tuning. 
In the next section we will show that upon adding further 
chains this line does in fact open into an extended region 
that develops into the Fibonacci phase in the 2D limit. 


IV. TOWARDS TWO DIMENSIONS 

Armed with our insights from the two-chain model ex¬ 
plored in the previous section, we wish to now deduce the 
phase diagram for our parafermion model in the full 2 D 
limit. For this purpose we use extensive DMRG simula¬ 


tions to complement earlier analytical studies that apply 
only in the case of weakly coupled chains (ts ^ ^1^2)- We 
are particularly interested in tracking the extent of the 
Fibonacci phase as the interchain coupling increases away 
from this limit, though we will attempt to address the 
nature of nearby states as well. To approach two dimen¬ 
sions we study the Hamiltonian of Eq. (§ on cylinders 
with Ny = 4 , 6, 8, and 10 sites around the circumference. 
The cylinders will be taken infinitely long in the hori¬ 
zontal direction, in part to avoid complications arising 
from possible gapless edge states. For our calculations 
we typically retain up to m = 5500 states in DMRG for 
truncation errors up to 10“^ and often as small as 10“^^. 


Fibonacci Phase on an Anisotropic Triangular 
Lattice 


Until specified ot herwis e we assum e ti^2,3 > 0 and fix 
and 


ti = t2 = t±_ (Secs. 


IV B 


IV C 


relax these assump¬ 


tions). Our goal here is to numerically assess how the sys¬ 
tem evolves as we vary to tune from the decoupled- 

chain limit = 0 up through the isotropic triangular 
lattice point t±_ = ts and beyond. From the coupled-wire 
analysis, one expects the 2 D system to immediately enter 
the gapped Fibonacci phase upon turning on any small 
but finite t±/ts 1 . In this limit the system should 
exhibit a quite long correlation length in the horizontal 
direction but—crucially—an arbitrarily short correlation 
length along y. This feature is extremely attractive for 
DMRG: it implies that infinite cylinders even with rel¬ 
atively small Ny possess local properties closely emulat¬ 
ing those of the fully 2D system of interest. Our sim¬ 
ulation results, even for cylinders as small as Ny = 4 , 
indeed show strong evidence that the system realizes the 
Fibonacci phase once > 0 . In what follows we will 
numerically recover the characteristics (i) through (1;) of 
the Fibonacci phase delineated near the end of Sec. [Hj 
For a wide range of t^/ts we observe two quasi¬ 
degenerate ground states—both in the triality 1 quantum 
number sector—obtained by starting the infinite DMRG 
algorithm in randomized initial states. (Which quasi¬ 
degenerate ground state appears depends on the pre¬ 
cise initial state used.) As we will argue later based on 
entanglement-entropy scaling, these two states span the 
entire ground-state subspace. Anticipating Fibonacci- 
phase physics, let us call the quasi-ground state with 
lower entanglement entropy | 1 ) and the other \s). Fig- 
urepU] shows the relative energy splitting — Ei)/\Ei\ 
of these states as a function of at various system 

sizes. For a fixed circumference Ny and sufficiently weak 
we observe a very small splitting (< 0.1%) with 
weak t± dependence. Beyond an V^-dependent scale of 
t_L, however, the splitting grows rapidly with the inter¬ 
chain coupling. We expect that the crossover between 
these behaviors transpires when the correlation length 
in the y direction becomes comparable to the circumfer¬ 
ence; this interpretation is consistent with the enhanced 








0.015 



FIG. 10. Relative energy splitting between quasi-degenerate 
ground states of cylinders with Ny = 4-10 as a function of 
t± (in units where ts = 1). Notice that for Ny > S the states 
remain nearly degenerate even up to the isotropic triangular 
lattice point t± = 1. This is the first evidence suggesting 
that in the 2D limit the Fibonacci phase enjoys a wide sta¬ 
bility window extending from weakly coupled chains past the 
isotropic triangular lattice limit. 


robustness of the degeneracy evident in Fig. upon in¬ 
creasing Ny. 

Next we examine the entanglement spectrum of both 
quasi-ground states, which for a topological phase is ex¬ 
pected to reveal the gapless low-energy spectrum in the 
presence of an open edge.^^ Following the approach of 
Ref. [ 62 l we plot the entanglement spectrum after apply¬ 
ing a non-universal overall shift and rescaling such that 
the lowest two levels match the spectrum of the chiral 
Z3 parafermion CFT. (Each entanglement “energy” is 
defined in terms of a reduced-density-matrix eigenvalue 
Pi through Ci = —\npi.) As shown in Fig. 
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for Ny = 4 

and = 0.2, all the remaining entanglement ener¬ 

gies exhibit the same pattern of degeneracies and relative 
energy-level spacings as the excitations of the chiral Z3 
parafermion CFT on a ring. Specifically, we find that 
the entanglement spectrum of state | 1 ) matches the CFT 
level pattern for the superselection sector corresponding 
to the primary fields l,'^,'^'*' (and their descendants), 
while the state \e) corresponds to the fields 5, a, . From 
the fusion algebra of the chiral fields, this suggests that 
1 5) and | 1 ) respectively carry Fibonacci and identity flux. 
We observed similar agreement of the entanglement spec¬ 
trum with the CFT prediction for other (small enough) 
values of t±/ts and for Ny = 6, 8,10, but choose to show 
the Ny = 4 results to emphasize how quickly the 2 D 
behavior sets in for these anisotropic cylinders. 


For further evidence that we are observing a topolog¬ 
ically ordered phase, we compute the topological entan¬ 
glement entropy 7^ for quasi-ground states \n = 1,^). 
For a “vertical” entanglement cut dividing our infinite 
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FIG. 11. Entanglement “energies” of the quasi-ground states 
|1) (left column) and \e) (right column) for the Ny = 4, 
t±/t 3 = 0.2 cylinder after shifting and rescaling the spectra to 
match the lowest two levels of the Z 3 parafermion GFT pre¬ 
diction. The numbers in parentheses indicate the predicted 
degeneracies. The same rescaling was used for both spectra, 
but with different shifts; thus the fit requires only three fitting 
parameters. 


cylinders into two semi-infinit e halv es, the entanglement 
entropy is predicted to scale 

Sn = aNy - 7n H- (8) 

where the first term encodes the boundary law expected 
to hold within a gapped phase. Because topologically or¬ 
dered states are locally indistiguishable in the 2 D limit, 
the coefficient a should be independent of the ground 
state n (although it is sensitive to other microscopic de¬ 
tails). Eurthermore, because DMRG is somewhat biased 
to low-entanglement states, it generally favors ground 
states of topological phases on infinite cylinders hav¬ 
ing well-defined topological flux through the cylinder, 
known as minimum entropy states . 1 ^ Eor such states, 
jn = ^og{V/dn), with dn the quantum dimension of 
the anyon type associated with the ground state 
and V = is the total quantum dimension of 

the theory. Notice that = 1 when summed 

over n, providing a way to check whether a complete set 
of ground states has been found. 

Eigure[T^fits our numerically computed entanglement 
entropies to the form in Eq. ^ . One subtlety here arises 
from the fact that the parafermions in our system can not 
appear in vacuum, but rather require a host system— 
i.e., a u = 2/3 quantum Hall state or similar ‘parent’ 
topological phase. Consequently, to back out the of 
interest we must shift our entanglement measurements by 
the topological entanglement entropy of the parent phase. 
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FIG. 12. Entanglement entropy fits for Ny =4,6,8 cylinders 
with two different interchain couplings t± (with ts = 1). All 
entanglement-entropy values are shifted downward by log 
as explained in the text. The y-intercepts closely match the 
topological entanglement entropies 7n predicted for the Fi¬ 
bonacci phase. 



71 

Is 

g-2-yi g- 27 . 

Exact 

logP 0.6430 

log(D/¥’) « 0.1617 

1 

t± = 0 . 4 ^ 

0.6235 

0.1393 

1.0442 

t± = 0.4*’ 

0.6306 

0.1538 

1.0186 

t± = 0.6 

0.6498 

0.1562 

1.0043 


TABLE I. Intercept values extracted from the fits in Fig. 
for t± = 0.4, 0.6 and ts = 1. Both agree with the theoretically 
predicted topological entanglement entropy shown in the ta¬ 
ble’s first row, and demonstrate that the set of ground states 
is complete. (The latter conclusion follows from the fact that 
the data in the right column are very close to unity.) For 
the t± = 0.4 system we give both (a) the intercepts shown 
in Fig, obtained from fitting all three Ny = 4,6,8 points 
and (b) intercepts obtained from only fitting the Ny = 4, 6 
entropies that could be computed more accurately. 


FIG. 13. Entanglement entropy difference between the quasi¬ 
degenerate ground states of Ny = 4,6,8,10 cylinders as a 
function of interchain coupling t± (with ts = 1). The hori¬ 
zontal dashed line denotes the thermodynamic-limit predic¬ 
tion \og(p ~ 0.481. Our Ny = 8 and 10 data further indicate 
that the Fibonacci phase survives even beyond the isotropic 
triangular lattice point in the 2D limit, corroborating the ev¬ 
idence presented in Fig. 


The numerical results for — Si in Fig. p! 3 | show good 
agreement with this value, especially for larger and/or 
more anistropic cylinders which are expected to have the 
weakest finite-size effects. 

Our numerical evidence regarding the entanglement 
entropy, entanglement spectra, and ground-state degen¬ 
eracy together strongly indicate the onset of a Fibonacci 
phase over a wide range of parameters.!^ These results 
not only corroborate analytical findings for the strongly 
anisotropic limit with <C 1; quite remarkably, 

Figs. [To] and [T^ also reveal that the Fibonacci phase per¬ 
sists into the isotropic-triangular-lattice case = 1 

and beyond! 


“7(2/3) = — log\/ 3 . (The unshifted 7^ values turn out 
to be negative.) Table [I] shows that after applying the 
shift, the observed 7^’s for two different magnitudes of 
t± agree very well with the theoretical prediction for the 
Fibonacci phase quoted in Sec. [TT| di = 1, = (/:? and 

V = where (/? = (! + a/ 5)/2 is the golden ratio. 

Furthermore, is very close to unity for 

both tj_ values, allowing us to deduce that the system 
admits no further ground states beyond |1) and j^). 

The difference — Si of entanglement entropies for 
the ground states has been found to converge more 
rapidly as a function of Ny than linear fits of each in¬ 
dividual Within the Fibonacci phase, we expect 

Ss - Si = -7^ +71 =\og(p ^ 0.481 ( 9 ) 

provided the cylinder size exceeds the correlation length. 


B. Extent of the Fibonacci Phase 

The previous subsection reported substantial evidence 
that the model in Eq. (§ realizes the Fibonacci phase 
along the line ti = ^2 = for a wide range of In 

light of this finding it is interesting to now explore the 
extent of the Fibonacci phase for general > 0 . 

To address this question we fix the ratio ti/ts and then 
vary ^2 from 0 to ti, thus scanning a ray in parameter 
space. Along this line we compute the bipartite entangle¬ 
ment entropy of infinite cylinders, observing clear peaks 
in the entropy as a function of ^2 (Fig. p! 4 | ) that indi¬ 
cate a transition out of the Fibonacci phase (smoothed 
into a crossover due to finite-size effects; we provide sup¬ 
porting evidence for this interpretation below). These 
peaks are quite broad for systems with ti/ts ~ 1 but 
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t2 


FIG. 14. Peaks in Ny = 8 cylinder entanglement entropies as 
a function of t 2 for various fixed ti (ts = 1 throughout). To 
estimate the peak locations, which in the 2D limit coincide 
with the transition out of the Fibonacci phase, data points 
near the peaks were fitted to a quadratic form. Fits are shown 
as thick gray curves while the resulting peak locations are 
indicated by vertical dashed lines. The peak locations give 
the four Ny = S data points in Fig. The ti = 0.4, 0.5, 0.6 
data were computed using DMRG with a fixed number of 
states m. The ti = 0.3 curve was found by extrapolating 
fixed t 2 entropies as a power law in m. 


sharpen considerably as ti/t^ approaches zero, presum¬ 
ably because for small the ^-direction correlation 

length remains below the cylinder circumference except 
for points very close to the phase boundary. 

To estimate the locations of the entropy peaks, data 
points near each peak were fit to a quadratic as shown 
in Fig. 14 for the case of Ny = 8 cylinders. Figure [T^a) 
shows the locations of the entropy peaks thus obtained 
for Ny = 4 ,6, 8 cylinders. For reference the black G = ^2 
curve represents the critical line for the two-chain limit 
discussed in Sec. |ml we now see that, as expected, this 
line broadens into an extended phase on larger cylinders. 
As a useful consistency check we note that for a 2D sys¬ 
tem composed of weakly coupled chains one can analyt¬ 
ically constrain the shape of the phase boundary, which 
should precisely coincide with these entropy peaks ex¬ 
trapolated to the Ny oo limit. In particular, scaling 
arguments reviewed in Appendix [P] predict that close to 
decoupled chains the critical couplings tic and t^c satisfy 


i2c -txc = C {t2c + 


( 10 ) 


where C denotes a non-universal constant and we have 
employed units where ts = 1 . Fits of the data to this 
form appear in Fig. [T^a) (solid lines). Despite having 
only one fitting parameter—and using data for moder¬ 
ately coupled chains—the fits are remarkably good even 
for the smaller cylinders. The agreement motivated us to 
also fit the data to an arbitrary power law; we find that 
for each Ny the fitted exponent agrees with the predicted 



t2 



FIG. 15. (a) Estimated extent of the Fibonacci phase near 
the weakly-coupled-chain limit. Data points correspond to 
the locations of peaks in the entanglement entropy observed 
when tuning 0 < t 2 < ti for fixed ti in Ny = 4, 6, 8 cylinders 
(ts = 1 throughout). The solid curves are fits of the data to 
the function (t 2 — ti) = C{t 2 -\-ti)^^^ with C an A^-dependent 
fitting parameter. Dashed curves show the fits reflected across 
the line ti = t 2 , since in the thermodynamic limit the Hamil¬ 
tonian is symmetric about this line, (b) Ny = 8 fit reflected 
under all ti-t 2 -tz permutations, estimating the full extent of 
the Fibonacci phase. The black squares and triangle mark the 
location of the isotropic square- and triangular-lattice models. 


value of 8/5 to within 10 %. 

In the 2 D limit, the phase diagram for Eq. (§ must 
be symmetric under permutations of ti, ^2, and We 
can exploit this symmetry to roughly estimate the full 
extent of the Fibonacci phase for all ti,t2,t3 > 0 based 
on the DMRG data discussed above. As a first step, 
the dashed curves in Fig. [T^a) reflect our DMRG data 
about the ti = ^2 line. We note, however, that finite- 
systems break the permutation symmetry; this naive re¬ 
flection thus represents an approximation whose validity 
increases with Ny. Ideally we could have directly mea- 
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FIG. 16. Global phase diagram (som e assembly required) of the Hamiltonian in Eq. 0 with couplings constrained to 
\ti\ + \t 2 \ + ItsI = 1- Note that Fig. M^b) shows the face with ti,t 2 G 3 > 0 containing the Fibonacci phase. The six corners 
of the octahedron correspond to various exactly solvable decoupled-chain limits, and directions (a)-(e) mark weakly-coupled 
chain limits where DMRG shows strong evidence of the semion phase. Direction (f) indicates anisotropic-square-lattice systems 
studied in Fig, and small squares mark locations of the isotropic square lattice model in the 2D limit. 

Bottom right corner: Phase diagram after assembly; refer to YouTube for instructions. 


sured entropy peaks near the expected phase transition 
curve on the t2 > h side also. But we found finite cylin¬ 
ders within this region to be much more highly entangled, 
preventing an accurate DMRG study. We then further 
permute our Ny = 8 data to arrive at the phase diagram 
shown in Fig. H^b), parametrized in terms of a vector 

t = {ti,t2,t3) ( 11 ) 

with ti + t2 + ts = 1 and all couplings non-negative. 
Notably, the isotropic triangular lattice point sits deep 
within the center of the phase—consistent with our re¬ 
sults from the previous subsection. 


C. Global Phase Diagram of the ti—12—G Model 

Our numerical results presented above indicate that 
the isotropic square lattice (points in Fig. [T^b) marked 
by open squares) realizes a state distinct from the 
Fibonacci phase. Pinpointing the precise nature of 
this state is, however, rather nontrivial. Approaching 


the square lattice from the weakly-coupled-chain limit 
0 < ti/ts ^ 1 with t2 = 0 does not provide immediate 
insight, as the continuum limit of the ti coupling yields 
competing relevant interactions whose effect is difficult to 
understand analytically. Similarly, applying DMRG to 
small, anisotropic square-lattice cylinders does not pro¬ 
duce any clear evidence of spontaneous symmetry break¬ 
ing or signatures of known topological phases. States 
appearing when the model contains negative couplings 
also remain mysterious at this point, since throughout 
this section we have so far focused exclusively on regimes 
with ti 2,3 > 0 . If we allow for arbitrary signs of ^1,2,3, 
Fig. [Tbib) actually contains only one out of eight faces 
in the full phase diagram for the Hamiltonian in Eq. (§! 
(For instance, a second face arises if we take G < 0 but 
G,3 > 0 and so on; see Fig. for an ‘unfolded’ sketch 
of the eight faces.) In whatfollows we aim to fill in 
these gaps and deduce the ground states appearing in 
the broader parameter space. 

Fortunately, we have not yet exhausted all of the ana¬ 
lytically soluble windows: with negative couplings a new 
tractable weakly-coupled-chain limit emerges that lends a 
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FIG. 17. (a) States of the U(l)6 chiral CFT, organized by charge n. The highest weight states have energies ^ and thus form 

2 

a parabola. Their descendents have energies ^ + m for integers m = 0,1, 2,... with multiplicities given by the partitions of 
the integers 1,1, 2, 3,... . Combining all the states for even n yields the edge spectra for the |ls) ground state, shown on the 
left axis. The states for odd n constitute the edge spectra for the \s) ground state, shown on the right axis, (b) Numerical 
entanglement spectra of the quasi-ground states |ls) (blue, left column) and \s) (red, right column) at t = (—0.1, —0.1, —1) for 
the Ny = S cylinder after shifting and rescaling the spectra to match the lowest two levels of the U(l)6 CFT prediction. We 
used the same rescaling for both spectra but different shifts. Numbers in parentheses denote the predicted degeneracies. The 
excellent agreement with analytics strongly supports the onset of the semion phase. 


great deal of insight into the problem. Consider the case 
ti = t2 = 0 with negative intrachain comling ts < 0 . Un¬ 
der a Fradkin-Kadanoff transformation^ each decoupled 
chain here maps to the self-dual point of the antiferro¬ 
magnetic three-state Potts model, which is described by 
a non-chiral U(l)g CFT with iF-matrix iF = (^ _g) and 
central charge c = 1 . 1 ^ We respectively denote the cor¬ 
responding left- and right-moving chiral boson fields for 
chain y by (friy) and (j)R{y)^ both 27 r-periodic in our nor¬ 
malization conventions. 

Just as with positive ts, we can once again bootstrap 
off of these decoupled critical chains to deduce the sys¬ 
tem’s response to small but finite ti^2- To do so we must 
expand the interchain couplings in Eq. (§ in terms of 
low-energy fields captured by the CFT. Reference |68| al¬ 
ready provided the desired mapping between lattice op¬ 
erators and continuum fields for the self-dual antiferro¬ 
magnetic Potts model; the result reads 

an/L ~ + • • • (12) 

where the ellipsis denotes subdominant contributions. 
Using this mapping, the interchain Hamiltonian 


{h + + H.c. 

(13) 

- t2)X{4>L,R{y),(f>L,R{y + r) • 

The first term possesses scaling dimension A = | at 
the decoupled-chain fixed point. In the second line X 
represents (numerous) subleading terms with dimension 
A = I; these can be enumerated explicitly but we refrain 
from doing so here. 

Except with fine-tuning, the 2D system’s properties 
are generically dominated by the more relevant (ti + ^2) 
term above. Its effect is relatively easy to understand for 
systems composed of a large but finite number of chains. 
When (ti + ^2) 7^ 0 , the first line in Hj_ directly couples 
the right-mover of chain y with the left-mover of ^ + 1. 
This hybridization gaps out the entire bulk but leaves 
‘unpaired’ chiral left- and right-moving U(l)6 fields at 
the first and last chains, respectively. The boundary CET 
thus admits chiral central charge c = — 1 and contains 
primary fields e^'^^ with scaling dimensions A = ^ and 

spins (0 is the gapless 27 r-periodic chiral boson 

field at the edge and n is an integer ‘charge’). 


takes the form 

y 

+ if 
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Note that these conclusions hold independent of the 
sign of (ti + ^2). 

Crucially, this boundary CFT does not describe the 
edge with the vacuum since the lattice parafermion op¬ 
erators comprising our system reside in the interior of a 
fractionalized quantum Hall fluid. Some care is therefore 
necessary when using the bulk-boundary correspondence 
to deduce properties of the gapped phase we have en¬ 
tered. All fields in the boundary CFT must correspond 
either to bulk quasiparticles originating from the cou¬ 
pled parafermions of interest or quasiparticles native to 
the surrounding quantum Hall fluid. To determine the 
former we must therefore ‘mod out’ by the latter. For 
concreteness, consider the first chain where 0 = 0 l(1 ) 
describes the gapless degress of freedom. Parafermion 
trilinears aLaLaL are bosonic and thus e^'^^ is local given 
our operator identification in Eq. ( 12 ). 

While single parafermion operators aL are clearly not 
bosonic, their nontrivial statistics originates from the 
fractionalized bulk quasiparticle of the parent quantum 
Hall fluid; in this sense aL ^ should also be consid¬ 
ered ‘local’. The fields for n G Z therefore corre¬ 

spond to the quasiparticle content of the parent phase. 
Such fields are all related by some ‘local’ operator and 
belong to the same superselection sector; likewise, the 
fields with odd n form another superselection sector. 
One can succinctly organize the field content—modulo 
the trivial boson e ^^"^—into the sets {1, and 

{gzc/), gSzc/), g 52 (/)} Modding out further by the ‘local’ field 
e^'^^ leaves only 1 and corresponding to the quasi¬ 
particles arising from our coupled parafermions. We con¬ 
clude that for weakly coupled chains with ts < 0, the 
(ti +^2) term in Eq. ( 13 ) drives the system into a topolog¬ 
ically ordered phase with a two-fold ground state degen¬ 
eracy on a torus or infinite cylinder. Easing the boundary 
field e^'^^ with itself yields e^'^^ x e^'^^ ^ e^^^ ^ a trivial bo¬ 
son, hence the nontrivial anyon has quantum dimension 
d = 1 . Moreover, since the spin of the boundary field 

is = i, this nontrivial bulk anyon is a semion.l^ 

We will therefore refer to this Abelian state as the 
“semion phase”. 

We confirm the presence of the semion phase by first 
applying infinite-cylinder DMRG to the Hamiltonian 
with t = (—0.1, —0.1, —1), which lies along the direction 
labeled (a) in Eig. 16 Starting from random initial states 
our simulations indeed recover two quasi-ground states 
that we denote jig) and \s). (The ‘s’ subscript on the 
former is used to distinguish from the state |1) defined 
for the Eibona cci pha se.) Their entanglement spectra, 
shown in Eig. 17 (b)[ a gree w ell with the U(l)6 GET 
counting shown in Eig. | 17 (a) We also find very simi¬ 
lar results (not shown) for weakly-coupled-chain systems 
with ti = t2 > 0 along direction (b) of Eig. 16 , and for 
systems along direction (e). Note that in the 2D limit 
directions (b) and (e) are related under permutation of 
t2 and ts. 

The less-relevant {ti — ^2) term in Eq. ( 13 ) plays an 



FIG. 18. Relative energy splitting between quasi-degenerate 
ground states of cylinders with Ny = 4-10 as a function of 
ti (in units where ts = 1 with t 2 = 0 throughout). The 
horizontal axis corresponds to direction (f) in Fig. 


important role only when (ti +^2) is tuned near zero. As 
mentioned above, this part of the Hamiltonian contains 
a sum of competing terms, each with scaling dimension 
I, making it difficult to determine analytically which 2 D 
phase it favors. Nevertheless, by turning again to nu¬ 
merics we find that along the line ti = —^2 close to the 

the 


16 


t = ( 0 , 0 , — 1 ) point [directions (c) and (d) in Fig. 
system still exhibits two quasi-ground states with entan¬ 
glement spectra consistent with the semion phase. This 
observation—along with the results for directions (a), 
(b), and (e)—leads us to believe that the semion phase 
occupies a large swath of parameter space around the 
symmetry-related points t = (0,0,—1), t = (0, —1,0), 
and t = (—1, 0, 0). 

Let us now return to the question of which phase 
emerges in the isotropic square lattice limit with posi¬ 
tive couplings, that is, t = (1,0,1) and its permutations. 
Directly approaching the t = ( 1 , 0 , 1 ) point from the 
decoupled-chain limit along direction (f) in Fig. 16 is un¬ 


fortunately challenging for DMRG. Here the ^-direction 
correlation length seems to grow very quickly as a func¬ 
tion of ti/ts relative to the slower growth encountered for 
the anisotropic triangular case upon varying As 

indirect evidence of the growth of y correlations. Fig. [ 18 ] 
shows the relative energy splitting of the two anisotropic- 
square lattice quasi-ground states as a function of ti 
(with ^2 = 0 and ts = 1). Gompared to Fig. 10 for 
the anisotropic triangular lattice, the anisotropic-square- 
lattice ground states separate more rapidly with ti and 
remain significantly split up to larger system sizes. Nev¬ 
ertheless, we consistently observe two quasi-ground states 
on anisotropic square lattice cylinders up through the 
isotropic point ti = ts = 1 (however, we cannot rule out 
the possibility that more than two exist). Gombined with 
our strong evidence that the semion phase appears along 
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directions (b) and (e) in Fig.[l^and persists in a broad re¬ 
gion around t = (0,0, —1) and symmetry-related points, 
we find it reasonable to conjecture that the isotropic 
square-lattice system also lies within the semion phase. 

Finally, we briefly comment on the symmetry-related 
regions labeled with question marks in Fig. 16 These 


parameter regimes represent weakly coupled chains de¬ 
scribed by Z3 parafermion CFT’s hybridized via negative 
interchain couplings. Consider the special case where 
ts = 1 with small ti = ^2 = t_L < 0 . Here the most 
relevant term generated by t± corresponds to the pertur¬ 
bation that drives the c = 4/5 parafermion CFT to the 
trier it ical Ising theory with c = 7 / 10 . 1 ^ We believe that 
the resulting critical state is unstable within our model, 
but one would need to assess the behavior of other less 
relevant terms to draw any conclusions. The situation 
should be contrasted to the case with > 0 where the 
leading perturbation opens a gap directly and yields the 
Fibonacci phase. Determining the fate of the system here 
is therefore expected to be especially delicate, so we leave 
an exploration of this subtle regime to future work. 


V. CONCLUSIONS 

Coupling arrays of extrinsic defects—e.g., those that 
bind parafermionic zero modes—provides a promising av¬ 
enue towards realizing exotic phases of matter, particu¬ 
larly given possible realizations in quantum-Hall archi¬ 
tectures. In an early work, Burrello et al.^^ constructed 
a model of coupled parafermions supporting an Abelian 
topologically ordered state that generalizes the toric 
code. Reference | 43 ] explored an alternative setup con¬ 
sisting of weakly coupled parafermion chains that were 
analytically shown to enter a ‘Fibonacci phase’ with non- 
Abelian Fibonacci anyons that possess universal braid 
statistics. This study established proof of principle that 
one could assemble hardware for a universal topological 
quantum computer using well-understood Abelian phases 
of matter. We stress, however, that the quasi-ID limit 
invoked there served purely as a theoretical crutch. Ex¬ 
ploring the broader stability of the Fibonacci phase away 
from this soluble regime is therefore important for assess¬ 
ing the utility of this line of attack. Here we employed 
extensive DMRG simulations, bolstered by complemen¬ 
tary analytical results, to address precisely this issue in a 
parafermion model that contains decoupled chains and 
the isotropic square- and triangular-lattices as special 
cases. 

Reassuringly, our numerics reveal a remarkably 
resilient Fibonacci phase—evidenced by the telltale 
ground-state degeneracy, entanglement entropy, and 
entanglement spectrum—extending well beyond the 
weakly-coupled-chain limit; recall Fig. [T^b). Elsewhere 
in the phase diagram an Abelian topologically ordered 
state with semion excitations appears. Although we have 
not exhaustively studied the model’s parameter space, 
our simulations revealed only these two phases (which. 


curiously, comprise the only two topological orders with 
a single nontrivial quasiparticle type). Note that we 
also performed simulations of square-lattice systems with 
first- and second-neighbor interchain couplings, though 
for brevity the results were not presented. This alter¬ 
native setup arises by augmenting Eig. with upward- 
sloping diagonal bonds. We found a clear Eibonacci 
phase here as well, which nicely emphasizes the robust¬ 
ness of the physics against changes in microscopic details. 

It is worth re-emphasizing the exceedingly useful com¬ 
plementarity of coupled-chain analytics and our DMRG 
simulations. The quasi-ID limit predicts sharp finger¬ 
prints of the Eibonacci phase that, because of the ultra- 
short transverse correlation lengths characteristic of this 
regime, could be verified numerically even in systems 
composed of as few as four chains. This feature al¬ 
lowed us to controllably track the extent of the Eibonacci 
phase in regimes where we lack analytical control—which 
of course was the primary goal of our simulations. We 
expect that analogous quasi-ID deformations should be 
useful for numerics in a variety of other contexts as well. 
One other aspect of the model is also worth stressing: it 
consists entirely of quadratic near-neighbor parafermion 
couplings yet captures highly exotic topological phases 
of matter; this is possible because even with onl y b ilin- 
ear terms the Hamiltonian is strongly interacting.^ The 
model’s simplicity is certainly a boon for the prospect of 
realizing Eibonacci anyons in the underlying quantum- 
Hall setups. 

Our results also yield indirect implications for systems 
built out of local bosonic degrees of freedom. In an inter¬ 
esting recent study, Barkeshli et al.l^ introduced general¬ 
ized Kitaev honeycomb models for which the constituent 
spin operators could be decomposed in terms of ‘slave 
parafermions’. Within that representation the authors 
identified a quasi-ID limit that also realizes a phase with 
Eibonacci anyons, among other quasiparticles. Refer¬ 
ence [ 73 ] conjectured the appearance of Eibonacci anyons 
in related local bosonic models (though Abelian topo¬ 
logical orders were found in a mean-field treatment). 
By applying the ‘inverse’ logic from these works one 
can back out local bosonic setups corresponding to the 
parafermion system that we analyzed. Doing so imme¬ 
diately allows one to construct bosonic lattice models 
supporting a phase with Eibonacci anyons—though the 
Hamiltonians involve somewhat elaborate multi-boson 
interactions. It would be interesting in future numerical 
work to explore the extent of that state when the model 
is deformed towards a simpler limit with only quadratic 
couplings. 

One particularly intriguiging wide-open question is 
whether the Eibonacci phase can appear in a reason¬ 
able spatially uniform Hamiltonian describing Abelian 
quantum Hall/superconductor heterostructures or bilay¬ 
ers (i.e., with the periodic modulations trapping the 
parafermion modes ‘smeared out’). If found, such setups 
might offer relatively simple experimental routes to Ei¬ 
bonacci anyons. The fact that we identified an isotropic 
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triangular lattice model supporting the Fibonacci phase 
is encouraging and may offer some clues as to how one 
can construct spatially uniform analogues. We also note 
that Barkeshli and VaezP^ proposed remarkably simple 
= I + I bilayers as uniform Fibonacci anyon hosts 
based on an analysis of the thin-torus limit. A concerted 
numerical push in this direction seems a very worthwhile 
complement to the lattice simulations performed here. 
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Appendix A: Mapping to Potts chain 


Here we very brieffy review the mapping of a ID 
Z3 pa raferm ion chain to the three-state quantum Potts 
model.l^^l^ Consider the Hamiltonian in Eq. Q, re¬ 
peated here for clarity, 




H.c. 


(AI) 


where again uj = We have chosen to express 

(somewhat arbitrarily) the Hamiltonian using ‘right’ 
parafermion operators that satisfy = I, a\^j = 0?^ ^^ 
and 


(A2) 

One could equally well employ deft’ parafermion opera¬ 
tors whose commutator instead reads 






(A 3 ) 


In either representation, the mapping to the Potts 
model proceeds upon nonlocally decomposing the 
parafermions via 




^L,2n—1 — An’ 

m<n 


aR, 2 n = WC 7 „ Tm (A 4 ) 

m<n 

aL,2n =OJ^Crn H A'(^- 5 ) 

m<n 


On the right sides we have introduced local bosonic Potts 
operators an , that satisfy 

= = 1, A = al, tI = tI (A6) 

together with the commutator 




(A 7 ) 


(In contrast to the parafermion operators, am and 
commute for m ^ n.) Notice that the right and left rep¬ 
resentations differ primarily in the form of their non-local 
Tm strings. One can r eadi ly verify that the decomposi¬ 
tions in Eq. ( |A 4 | ) and (A 5 ) preserve the parafermion op¬ 
erator algebra. Inserting these expressions into Eq. 
yields the desired quantum Potts Hamiltonian, 


H =-'^\j (A+lO'n + H.C.) + h{Tn + tI) 

n 


(AS) 


Appendix B: DMRG path 

When working with a model of coupled parafermions 
such as Eq. (§, it is necessary to specify the parafermion 
operator commutation relations by ordering the sites of 
the system. Such an ordering also specifies how the 
parafermion operators map into clock-model operators 
under a Eradkin-Kadanoff transformation. 

As an example, for a 2D array of parafermion operators 
a common choice of ordering is for operators on chain y 
to precede those on chain y' > while on the same chain 
operators remain ordered from left to right as in the ID 
case. Eor parafermion operators on the same chain y this 
implies 


aR,j{y)(^R,y{y) = e*( 2 -/ 3 ),(y) 

(Bl) 


while for operators on different chains, 

aR,j{y)aR,j’{y') = aRj,{y')aRj{y) 

aL,jiy)aL,yiy') = aL,y{y')aL,j{y) 

(B 2 ) 

regardless of the order of j and j' . When simulating a 2 D 
system using DMRG, which requires choosing a ID path 
for the matrix-product-state wavefunction, it is conve¬ 
nient to choose the DMRG path to follow the parafermion 
ordering. However, for DMRG to follow the above order¬ 
ing would lead to a very steep growth in computational 
cost for systems with a significant number of sites along 
the X direction. Eortunately, however, other physically 
equivalent parafermion orderings are possible. 

The alternative ordering (and corresponding DMRG 
path) shown in Eig.[^is reasonably efficient for perform¬ 
ing DMRG calculations. It also makes implementing a 
Hamiltonian such as Eq. relatively simple since none 
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FIG. 19. Possible ordering of parafermion operators in an 
Ny = 4 cylinder. For the Hamiltonian in Eq. ([^, this choice 
yields no interactions that cross the parafermion ordering 
path. Parafermion ‘hopping’ processes that pass through re¬ 
gions labeled Rot L are expressed in terms of the correspond¬ 
ing operator representation, e.g., or a^j^agL- 



FIG. 20. Ordering of parafermion operators used in DMRG 
for the present work (shown for a cylinder with Ny = 4). Note 
that in terms of the underlying clock degrees of freedom seen 
by DMRG, this path requires a circumference of only two sites 
in the y direction. 


of the interactions cross the parafermion path. To map 
an interaction involving, say, sites 3 an d 6 into clock- 
operator form, one uses the mapping Eq. (A4), using the 
R representation as this process ‘hops’ parafermions un¬ 
der the path (through the region labeled R in Fig. 19). 
The resulting correspondence is 



FIG. 21. ‘Braiding’ parafermion sites 3 and 4 clockwise two 
times transforms the pairwise interaction between sites 3 and 
7, defined originally in the L representation, into an interac¬ 
tion that crosses the parafermion ordering path. This scheme 
allows us to simulate the model of Eq. § using the numeri¬ 
cally very convenient ordering shown in Eig. |20[ 


sites along the ^-direction by a factor of two compared 
with a path that pairs parafermions into clock sites along 
the x-direction. 

Relative to the conventions specified in Fig. im¬ 
plementing the model of Eq. <§ using the ordering in 
Fig. is subtler since here many Hamiltonian terms 
involve processes that cross the parafermion path (e.g., 
coupling between sites 3 and 7). The clock-operator form 
of such processes can be derived from a simpler ‘reference’ 
parafermion ordering in which the pairwise interaction 
involving the same two sites does not cross the ordering 
path. One way to carry out this derivation is to apply a 
unitary transformation that performs a clockwise ‘braid’ 
of the parafermion sites (j, j + 1) as follows 


= CrjT2T3Cr3 . (B3) 

The factor of uj is included above so that the Hamiltonian 
term involving this bilinear +H.c.) is symmet¬ 

ric under charge conjugation, defined as the mapping: 

O^R/L ^ 

However, because of the way short-range interactions 
in 2D map to longer-ranged interactions in the ID path 
used by DMRG, the most efficient path choice for DMRG 
corresponds to ordering the sites by columns such that 
all operators in column x are ordered before those in 
X + 1. Therefore in the present work we we have used 
the ordering illustrated in Fig. Because such a path 
‘pairs’ parafermions into clock sites along the ^-direction, 
it has the added benefit of reducing the number of clock 


^ 2 t t (B5) 

+ l ^ ^ ^R,j + l^R,j 

for parafermionic operators in the R representation; for 
those in the L representation 

The above mapping corresponds to the braid transfor¬ 
mation derived for parafermionic zero modes in Ref. |32l 
As an illustrative example, consider the interaction 
‘hopping’ a parafermion from site 7 to site 3 in Fig. 
which crosses the ordering path. We can construct such 
an interaction beginning from the simpler process de¬ 
picted in Fig. [2^a) that hops a parafermion from site 
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FIG. 22. Parafermion ordering path used to discuss the self¬ 
duality of the two-leg ladder model studied in Section |III| 
Solid lines and numbering of the parafermion operators in¬ 
dicate the ordering path, with coincides with the ti and t 2 
interactions. The ts interactions, shown as dashed lines, in¬ 
volve processes that take place above or below the path. 

7 to site 3 while staying above the path, thus involving 
parafermion operators in the L representation: 


Appendix we show that the ladder exhibits a variant of 
self-duality when ti = ^2, independent of the strength 
of ts, that sharply constrains the location of the critical 
line separating the topological and trivial phases found 
numerically; recall Fig.[^ Following the discussion of Ap¬ 
pendix it will prove useful to equivalently rewrite the 
Hamiltonian in terms of a single chain of parafermions 
ordered along a path that follows the ti and t 2 bonds. 
The specific path, including conventions for R and L, 
appears in Fig. 

In this convention the two-leg ladder Hamiltonian be¬ 
comes 

3 

+ + H.C.] . 

(Cl) 


uja^^aLj = . (B7) 

(We display this operator in both parafermionic and 
clock-operator form for clarity.) Braiding sites 3 and 4 
clockwise maps ^ transforming this interac¬ 

tion into the one depicted in Fig. [^b) as follows: 

cjcrt nOiL 7 —^ cjtrt aC^l 7 
L,3 L ,4 R 7 

= cr^T^a4 . 

Finally, braiding sites 3 and 4 clockwise once more maps 
o^L ,4 3, yielding 


The first line coincides with the nearest-neighbor 
parafermion chain discussed in Sec. [HI and Appendix [A] 
the second line simply adds second-neighbor interactions. 

Equations (A4) and (A5) allow us to alternatively ex¬ 
press the Hamiltonian in Potts language: 


ff2-leg = - + H.C.) + ti{Tj + t]) 

3 

+ + H.c.]}. 


(C2) 


Next we introduce a duality transformation, defining op¬ 
erators 


oja\^^^aLj -)> aLj 

= <^L,3<^L,4<^L,7 (B9) 

2 t t 

= UJ (J2T2T^CF4 


The key result here is the last line above, w hich gives 
the clock-operator form of the interaction Eq. (|B7[ ), but 
consistent with the ordering convention of Eig. [2Q| 
Because the Eig. path is a less natural, but nu¬ 
merically very convenient, choice for the Hamiltonian 
Eq. (§ in that interactions cross the parafermion path, 
we confirmed that our mapping to clock operators passed 
a number of consistency checks. One of these is that 
parafermion bilinears which do not intersect should com¬ 
mute. Another is that the resulting clock-operator 
forms of each interaction term should respect charge- 
conjugation symmetry. We also checked that imple¬ 
menting the Hamiltonian Eq. via both the Eig. 


and Eig. paths gives the same ground-state energy in 
DMRG for a wide range of system sizes and Hamiltonian 
parameters. 


Appendix C: Self-duality in the Two-leg Ladder 

Consider the model of Eq. <§ defined on a two-leg lad¬ 
der (i.e., Ny = 2) with ti^2,3 purely real as usual. In this 


H= d-iv 

i<j 

that satisfy the same algebra as the original a, r variables, 
e.g., RjUj = ujUjRj. In this dual representation we obtain 


^2-leg — - + ^2(z7j + Z^j) 

3 

+ i3[w(d + + H.c.]}. 


(C4) 


The first line retains the same form as in Eq. (C2) but 


with ti and t 2 swapped. Crucially, the ts term in the 
second line also retains the same form upon additionally 
performing an antiunitary ‘time-reversal’ transformation 
that sends 


Uj , fJjj jJjj. (C^) 

(Note that this transformation preserves the commuta¬ 
tion relations for the dual variables.) Thus, as adver¬ 
tised, i72-ieg is invariant under duality followed by time 
reversal when ti = t 2 for arbitrary This observa¬ 
tion constrains the critical line separating the trivial and 
topological gapped phases for the two-leg ladder to reside 
precisely along the line ti = ^2, which is indeed borne out 
in our numerical simulations (Eig. |^. 
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Appendix D: Scaling argument for Fibonacci phase 
boundary 


With ti = t 2 = 0 and ts > 0 our Hamiltonian in 
Eq. describes Ny independent parafermion chains, 
each tuned to the self-dual critical point described by a 
non-chiral Z3 parafermion CFT. This CFT contains six 
primary fields {1,7/;,'0'*', 5, a, cr’*'} with respective scaling 
dimensions {0, |, Reference [53l relates these 

fields directly to the lattice parafermion operators of our 
model in the long wavelength limit: 


o^Rj ^ a^jjR + b{—iy^cjReL + • • •: 
aL,j ^ a^pL + + ..., 


(Dl) 


with non-universal coefficients a, b and ^ab denoting the 
fusion product of fields A and B. The first term in each 
expansion has a scaling dimension of |, while the second 
has a (smaller) scaling dimension of The ellipses in¬ 
dicate higher-order terms with scaling dimension larger 
than 1. 


Using Eqs. (Dl), the interchain Hamiltonian of Eq. ^ 
becomes 


^{ti+t2)Loa^ipiij,y)^Rij,y + l) 

3 

+ (il 1) +H.C. 

+ ... (D2) 


where we have retained only the two most relevant cou¬ 
plings. The wire construction of the Fibonacci phase in 


which we are interested^ arises from the (ti + ^2) term 
in the above expansion. When ti = t 2 > 0, this cou¬ 
pling clearly dominates. In this case, the left movers in 
the CFT of one chain couple only to the right movers of 
the next, gapping out the interior modes and leaving a 
chiral Z3 parafermion CFT at the top and bottom edges. 
However, it is important to note that the second term 
(for which left and right movers of each chain interact, 
leaving the fate of the system less clear) is actually more 
relevant: The first coupling has scaling dimension 4/3, 
while the second has dimension 14/15. 

One can use this scaling information to estimate the 
shape of the phase boundary separating the Fibonacci 
phase from the state favored by the (ti — ^2) term. Defin¬ 
ing X± = ti±t2, the renormalization group flow equations 
for the couplings A± read 

diX± = (2 - A±)A±, (D3) 

where I is a logarithmic rescaling factor while A+ = | 
and A_ = ^ denote the corresponding scaling dimen¬ 
sions. By rearranging the flow equation, one can see that 
the combination does not change with 

the renormalization flow. At criticality the A+ and A_ 
terms compete to a draw. The phase boundary should 
therefore reside along a constant value of this flow invari¬ 
ant in order for the two coupling constants to remain in 
competition as I increases. In our case, this means that 
the phase boundary should occur approximately along a 
curve given by 


he -hc = C {he + hef/^ (D4) 

for some non-universal constant C. 
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